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Oh' Abstract 

|' With present and upcoming SUSY searches both directly, indirectly and at acceler- 

ators, the need for accurate calculations is large. We will here go through some of the 
tools available both from a dark matter point of view and at accelerators. For natural 
reasons, we will focus on public tools, even though there are some rather sophisticated 
private tools as well. 

> 

00 . 

t^J- , The long awaited Large Hadron Collider (LHC) is expected to start taking data in 2009. 

■ The LHC research program has traditionally been centered around the discovery of the Higgs 

boson. However, the standard model description of this particle calls for New Physics. Until a 
£C) • few years ago, the epitome of this New Physics has been supersymmetry, which when endowed 

with a discrete symmetry called R-parity furnishes a good dark matter candidate. Recently 
a few alternatives have been put forward. Originally, they were confined to solving the Higgs 
problem, but it has been discovered that, generically, their most viable implementation (in 
accord with electroweak precision data, proton decay, etc.) fares far better if a discrete 
symmetry is embedded in the model. The discrete symmetry is behind the existence of a 
possible dark matter candidate. 

From another viewpoint, stressed in many parts of this book, the last few years have 
witnessed spectacular advances in cosmology and astrophysics confirming that ordinary mat- 
ter is a minute part of what constitutes the Universe at large. At the same time in which 
the LHC will be gathering data, a host of non collider astrophysical and cosmological ob- 
servations with ever increasing accuracy will be carried out in search of dark matter. For 
example, the upcoming PLANCK experiment will make cosmology enter the era of precision 
measurements akin of what we witnessed with the LEP experiments. 

The emergence of this new paradigm means it is of utmost importance to analyse and 
combine data from these upcoming observations with those at the LHC. This will also pave 
the way to search strategies for the next Linear Collider, LC. This important program is only 
possible if a cross-border particle-astroparticle collaboration is set up having at its disposal 
common or complementary tools to conduct global searches and analyses. Many groups, from 
erstwhile diverse communities, are now developing, improving, generalising, interfacing and 
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exploiting such tools for the prediction and analysis of dark matter signals from a combination 
of terrestrial and non terrestrial observations, paying particular attention to the astrophysical 
uncertainties. Most of this work has been conducted in the context of supersymmetry, but 
the latest numerical tools are not limited to it. 

In this chapter, we will go through some of the tools available both from a dark matter 
point of view and at accelerators. For natural reasons, we will focus on public tools, even 
though there are some rather sophisticated private tools as well. For supersymmetric dark 
matter calculations, one of the first public tools available was Neutdriver [HE]- It was a 
precursor in the field, but has by now been superceded by other more sophisticated tools. 
There are currently three publicly available codes for calculations of dark matter densities 
and dark matter signals: DarkSUSY 010, micrOMEGAs [BJ [71 HI H] and IsaRED (part of 
ISASUSY/Isajet [10]) 

1 Annihilation cross section and the relic density 

The general theory behind relic density calculations of dark matter particles is given in 
Chapter 7 . Here we will focus on supersymmetric dark matter (neutralinos) and various 
tools available for calculating the relic density. There are currently three publicly available 
codes for calculating the relic density of neutralinos: DarkSUSY, micrOMEGAs and IsaREld. 
All three of these codes are capable of reading (and sometimes writing) SUSY Les Houches 
Accord (SLHA) files [T31II1] which allows for an easy interface between these codes and other 
tools to be described in Section 

We will in the following refer to these codes and how they calculate the relic density. We 
will use the notation of Chapter 7 [11) and only write down the equations needed to facilitate 
our discussions here. 

1.1 The Boltzmann equation 

In most supersymmetric models of interest for dark matter phenomenology, the lightest 
neutralino, \ is the lightest supersymmetric particle and our dark matter candidate. As 
such, we want to calculate the relic density of neutralinos in the Universe today as accurately 
as possible, which means that we need to solve the Boltzmann equation. 



where n x is the number density of neutralinos, H is the Hubble parameter, (a ann v) is the 
thermally averaged annihilation cross section and n Xt oq is the equilibrium number density of 
neutralinos. This equation needs to be solved over time (or temperature) properly calculating 
the thermal average at each time step. When the neutralinos no longer can follow the chemical 
equilibrium density n XtCq , they are said to freeze-out. There are several complications in 
solving Eq. (fT|); for example, we may have resonances and thresholds in our annihilation cross 
section. The solution to this is to calculate the annihilation cross section in general relativistic 
form, for arbitrary relative velocities, v. Another complication is that other supersymmetric 
particles of similar mass will be present during freeze-out of the neutralinos. To solve this we 
need to take into account the so-called coannihilations between all the SUSY particles that 
are almost degenerate in mass with the neutralino (in practice, it is often enough to consider 
coannihilations between all SUSY particles up to about 50% heavier than the neutralino). 

2 After this chapter was completed, a fourth code, Superlso Relic 1 1 21 has appeared as well. 
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Following the discussion in Chapter 7 [JTJ , we can solve for the total number density of SUSY 
particles, 



N 

(2) 



instead of only the neutralino number density. It is also advantageous to rephrase the Boltz- 
mann equation in terms of the abundance Y — n/s and use x = m x /T as independent 
variable instead of time or temperature. When coannihilations are included, the Boltzmann 
equation (JTJ) can then be written as 



dx V nMp 



1/2 1/2 

a eS v) (Y 2 - Yl) . (3) 



dY ( 45 \ 9* ™x , 



where (a c fiv) is given by [TTJ Eq. (7.18)] 



1.2 Solving the Boltzmann equation 

To solve the Boltzmann equation (j3j) we need to calculate the thermally averaged annihilation 
cross section (a e gv) for each given time (temperature). This is typically quite CPU- intensive, 
and we therefore need to use some tricks. In DarkSUSY [5], the solution is speeded up by 
tabulating W c s in [11] Eq. (7.19)] , but using the momentum of the x, p e g as independent 
variable instead of s. This tabulation takes extra care of thresholds and resonances making 
sure that they are tabulated properly. This tabulated W e g is then used to calculate the 
thermal average (<T e ftV) for each time (temperature), using [TTJ Eq. (7.18)] . The advantage 
with this method is that W e g does not depend on temperature and instead the temperature 
dependence of (a c gv) is completely taken care of by the other factors in [TTJ Eq. (7.18)]. 
Numerically, one needs to take special care of the modified Bessel functions K\ and Ki 
which both contain exponentials that need to be handled separately to avoid numerical 
underflows. The Boltzmann Equation ([3]) is then solved with a special implicit method with 
adaptive stepsize control, which is needed because the equation is stiff and develops numerical 
instabilities unless an implicit method is used. The details of the DarkSUSY method are as 
follows (see DarkSUSY manual 5 ). The derivative dY/dx in Eq. ([3]) is replaced with a finite 
difference (Y n +i—Y n )/h, where Y n = Y(x n ) and £„+i = x n + h. Then Y n+ i is computed 
in two ways: first, the right hand side of Eq. ([3]) is approximated with — ^[X n (Y 2 — Y e 2 „) + 

A„ + i(Y,j +1 — Y e 2 „ +1 )], where \{x) — (Ah/nMp) ^ 2 (gt^m/x 2 ) (a c gv), and an analytic 
solution is used for the resulting second-degree algebraic equation in Y n+ i; second, the right 
hand side of Eq. ([3]) is approximated with — X n +i(Y 2 +1 — Y 2 q „ +1 ) and an analytic solution 
of the algebraic equation for Y n +i is used. The stepsize h is reduced or increased to maintain 
the difference between the two approximate values of Y n +i within a specified error. Overall, 
the solution of the Boltzmann equation and the tabulation of W e s solves for the relic density 
to within about 1%. If needed, higher accurary can also be chosen as an option. 

In micrOMEGAs 7], (cr ann v) at a given temperature T is arrived at by performing a direct 
integration and does not therefore rely on a tabulation of the matrix elements squared. Two 
modes are provided to perform the integration. In the accurate mode the program evaluates 
all integrals by means of an adaptive Simpson integration routine. It automatically detects 
all singularities of the integrands and checks the precision of the calculation increasing the 
number of points until an accuracy of 10~ 3 is reached. In the default mode (fast mode) 
the accuracy is not checked but a set of points is sampled according to the behaviour of 
the integrand: poles, thresholds and Boltzmann suppression at high momentum. The first 
integral over scattering angles is performed by means of a 5-point Gauss formula. The 
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Figure 1: Total differential annihilation rate per unit volume dA/dp e s for the same model 
as in pTJ Fig. 7.3], evaluated at a temperature T = m x /20, typical of freeze-out. Notice the 
Boltzmann suppression at high p c g . 

accuracy of this mode is generally about 1%. The user can also test the precision of the 
approximation based on expanding the cross section in terms of its s and p-wave components. 

1/2 

In the Boltzmann Equation, we need to know g c ^{T) and h c s(T) that enter gj through 
[TT1 Eq. (7.9)]. In DarkSUSY, the default is to use the estimates in Ref. [TS], but other options 
are also available. Typically, different estimates of g e g(T) and h e g(T) translate into relic 
densities different by a few percent. 

IsaRED on the other hand does not solve the Boltzmann equation numerically, instead 
it finds the freeze-out temperature (the temperature when the annihilation rate equals the 
expansion rate of the Universe) and calculates the relic density from that (including remnant 
annihilations at later times). For the thermal averaged annihilation cross section, it uses the 
same relativistic treatment as DarkSUSY and micrOMEGAs. 

1.3 Coannihilation criteria 

In principle, one should include all SUSY particle coannihilations when calculating the relic 
density. However, the heavier they are, the less abundant they will be and can thus be 
neglected. This is important to speed up the calculation, as we would otherwise spend 
most of our CPU cycles on calculating non-important coannihilation cross sections. One 
can estimate which particles to include, by investigating (TTJ Eq. (7.18)]. The modified 
Bessel function K\ contains an exponential, the so called Boltzmann suppression, that will 
suppress all heavier particles. In Fig. [TJ we show dA/dp c g for the same model as in [TT1 
Fig. 7.3]. dA/dp e ff is essentially (apart from normalization) the integrand in the numerator 
in Eq. (7.18)]. Comparing the two figures, we clearly see the Boltzmann suppression 
of larger p e gf, i.e. heavier coannihilating particles. We can quantify this by comparing the 
Boltzmann suppression factor for two coannihilating particles with masses m, and rrij with 
the corresponding factor for the LSP, the x- The suppression factor for the coannihilating 
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particles compared to the x is roughly (neglecting the pl s in [TTJ Eq. (7.18)]) given by 



B = 



K 1 {{m l +m J )/T) 
K x {2m x /T) 



X 



X 



X 




(4) 



At freeze-out we typically have x ~ 20, which gives a suppression factor of B ~ 1CP 6 for 
coannihilation particles about 50% heavier than the LSP, \- in DarkSUSY, one can set the 
maximum mass fraction of coannihilation particles, / = m,i/m x that will be included in the 
calculation, whereas in micrOMEGAs one sets the minimum B to allow instead. The defaults 
in DarkSUSY (/ = 1.5) and micrOMEGAs (i? m i n = 10~ 6 ) are roughly equivalent. One should 
remember, though, that the value to choose depends on the particle physics model. For 
example, for chargino coannihilations, the x + X~ coannihilation cross section can be orders 
of magnitudes larger than the \X annihilation cross sections and one should choose / or B m j n 
so that one does not accidentally neglect coannihilations that are important. For the MSSM, 
the default values of DarkSUSY and micrOMEGAs are typically sufficient for all interesting 
cases. IsaRED instead includes a preset collection of particles that are of relevance for the 
mSUGRA setup. 

1.4 Annihilation cross section 

At the heart of the relic density calculation are the annihilation and coannihilation cross 
sections. In the MSSM there are over 2800 sub-processes (not counting charged-conjugate 
final states) that can in principle contribute in the relic density calculation. It appears at 
first sight to be a daunting task to provide such a general code. 

In DarkSUSY, all annihilation and coannihilation cross sections for the MSS]V0 are cal- 
culated at tree level by hand with the help of symbolic programs like Reduce, Form or 
Mathematica. The calculations are performed with general expressions for the vertices in the 
Feynman rules and the results are converted to Fortran code. The vertices are then calcu- 
lated numerically for any given MSSM model. The analytically calculated cross sections are 
differential in the angle of the outgoing particles, and the integration over the outgoing angle 
is performed numerically. 

In micrOMEGAs, on the other hand, any annihilation and coannihilation cross sections are 
calculated automatically and generated on the fly. This is possible thanks to an interface to 
CalcHEP [16], which is an automatic matrix element /cross sections generator. This automa- 
tion is carried one step further in that CalcHEP itself reads its MSSM model file (Feynman 
rules) from LanHEP [17] . which outputs the complete set of Feynman rules from a simple cod- 
ing of the Lagrangian, see section [5] In the first call to micrOMEGAs only those subprocesses 
needed for the given set of the MSSM parameters are generated. The corresponding "shared" 
library is stored on the user disk space and is accessible for all subsequent calls, thus each 
process is generated and compiled only once. This library is then filled with more and more 
processes whenever the user needs new processes for different MSSM scenarios. This avoids 
having to distribute a huge code with all the possible 2800 processes. 

Both methods have advantages and disadvantages. In the DarkSUSY setup, no recalcula- 
tion of the (analytical) annihilation cross sections is needed, which can speed things up. Also, 
the analytically calculated annihilation cross sections can be optimized to be faster. On the 
other hand, the micrOMEGAs setup makes it easier to adapt the code to non-MSSM cases. In 
both codes though, the actual Boltzmann equation solver is very general and works for any 
kind of WIMP dark matter, not only SUSY dark matter. In IsaRED, CompHEP [18] is used 
to calculate the annihilation and coannihilation cross sections for a subset of SUSY particles 
of relevance mostly for mSUGRA (the two lightest neutralinos, the lightest chargino, the 

3 Gluino coannihilations are currently not included. 
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left-handed eigenstates of sleptons and squarks, and gluinos). The expressions for the anni- 
hilation cross sections in IsaRED are not calculated on the fly, but are instead precalculatcd 
and included with the code. 



2 Direct detection 



Detailed expressions for detection rates in direct detection experiments are presented in 
Chapter 17 [19]. Here we focus on characteristics of the elastic scattering cross sections and 
event rates as they are implemented in numerical tools. 

Direct detection rates depend on the differential elastic WIMP-nucleus cross section 
d&WN I dEji , where Er is the energy of the recoiling nucleus: 



dE R . 
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(5) 



Here rriN is the nucleus mass, /in = m x mN j (m x + m^) is the WIMP-nucleus reduced mass 
(m x being the WIMP mass), v is the WIMP-nucleus relative velocity before the collision, 
Q.SI.SD are ^ e spin-independent and spin-dependent cross sections at zero momentum trans- 
fer, and Fgj sd (Er) are the squares of the corresponding form factors (also called structure 
functions). In terms of these quantities, the directional and non-directional direct detection 
rates read 
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(8) 
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v m in = [{ m NER)/{2y? N )] 1 ^ 2 is the minimum WIMP speed which can cause a recoil of energy 
Eft, fl q is the direction of the nucleus recoil momentum of magnitude q = v / 2mJv^R, and 
f(v m i n , fl q ) is the Radon transform of the velocity distribution function /(v). 

An important property of Eqs. © and (O is the factorization of the particle physics 
properties, uwn(Eh), and the astrophysics properties, pof(v m i n , Qq) and por](v m i n ). 

Dark matter codes such as DarkSUSY, IsaRED/RES and micrOMEGAs compute the particle 
physics and the astrophysics factors to various levels of precision and offer a number of choices 
for the form factors and the velocity distribution (more and more as they are upgraded) . All 
codes provide the zero-momentum transfer cross sections (although some are still limited to 
the axial and scalar couplings of supersymmetric neutralinos) , and most of the codes provide 
routines for direct detection rates off composite targets besides single nuclei. 

For the spin-independent part, dark matter codes use the factorized form a^ N (En) = 
Cq 7 Fgj(E R ) with, in the notation of Chapter 17 [T5] , 
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Various expressions for the spin- independent form factor Fgj(E R ) are typically available. 
For example, DarkSUSY 5 automatically selects the best available form factor among Sums- 
of-Gaussians, Fourier-Bessel, and Helm parametrizations (see [20] for a comparison of these 
approximations) . 

The spin-dependent part is often not factorized, so as to use the same functions provided 
by detailed simulations at zero and non-zero momentum transfer. With the by-now-standard 
normalization of the spin structure functions in |21j , one has 

a* D F§ D (E R ) = ^3^1 [ a ls 0Q (E R ) + alSn(E R ) + a oai S Q1 (E R )] (11) 

= ^ff^ KS PP (E R ) + a 2 n S nn (E R ) + a p a n S pn (E R )} , (12) 

where ciq = a p + a n , a\ — a p — a n , 

S PP (E R ) = Sqo(E r ) + Sn(E R ) + Soi(E R ) , (13) 
S nn (E R ) = S 00 (E R ) + S n (E R ) - S 01 (E R ) , (14) 
S pn (E R ) = 2[S 00 (E R )-S U (E R )}. (15) 

When the nuclear spin is approximated by the spin of the odd nucleon only, one finds [22) 
S PP (E R ) = — — — -, S nn (E R ) = 0, S pn (E R ) — 0, (16) 



for a proton-odd nucleus, and 



A^J(J + 1)(2J+1) 



S PP (E R ) = 0, S nn (E R ) = , S pn {E R ) — 0, (17) 

for a neutron-odd nucleus. Here Ajv is conventionally defined through the relation (iV|S|iV} = 
Atv(^V|J|^V), where \N) is the nuclear state, S is the nuclear spin, J is the nuclear total angular 
momentum. Tables of \ 2 N J(J + 1) values for several nuclei can be found in [23] and [2~I] . 

The quantities / p , f n , B Nl a , a x are sums of products of the WIMP-quark and WIMP- 
gluon coupling constants ctq'Q' A ' P ' T (for scalar, vector, axial, pseudoscalar, and tensor cur- 
rents) and of the contributions /tg, frq and A 9 of the gluons and each quark flavor to the 
mass and spin of protons and neutrons. Values for the nucleonic matrix elements of glu- 
ons and quarks, in practice values for /tg, frq and A 9 , are either hardcoded or settable 
by the user. Values for the effective coupling constants a s,v,A,p,T are gj^gj. precomputed 
analytically (DarkSUSY) or computed numerically on the fly (micrOMEGAs). For example, the 
effective lagrangian at the zero momentum transfer for the interaction of a fermionic WIMP 
X with quarks q reads 

£f = a s q xx qq + o% xi^x qi^q 

+ a q X757mX qibY'q + \a? XV^X qcr^q. (18) 

In the case of a Majorana WIMP, like the neutralino in the MSSM, only operators even 
under x ^ X are possible (i.e. = a? = 0). In micrOMEGAs, the numerical values of 
the coefficients a q are obtained combining appropriate matrix elements for xq ~^ Xq an d 
Xq ~^ Xq scattering at zero momentum transfer. For example, 

S, y ■ (q(Pi),x(P2)\SOs\q(pi),x(P2)) 
" 9 9 " 1 (q(pi),x(P2)\O s O s \ q ( Pl ), X (p 2 )) 



s_ v = ■ (q(Pi),x(P2)\SOs\q(pi),x(P2)) nq , 
q q (q(pi),x(P2)\OsO s \q( Pl ),x(P2)} 
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where Os = XX QQy the S'-matrix S — 1 — iC is obtained from the complete Lagrangian at the 
quark level, and the scattering matrix elements on the right hand sides are computed with 
CalcHEP. More general cases, including a generic local WIMP-quark operators and WIMPs 
with spin-0 and spin-1, are presented in [9]. 

Loop contributions are essential in the treatment of the WIMP-quark and especially 
WIMP-gluon coupling constants. For example, for neutral WIMPs like the supersymmetric 
neutralino, there is no neutralino-gluon coupling at the tree-level and the gluon contribution 
to a s arises at the one-loop level. Complete analytic one-loop calculations for neutralino- 
quark and neutralino-gluon couplings were performed in |25i I26j ; these formulas are incor- 
porated in DarkSUSY. Automatic numerical calculations of all a^' V ' A ' F ' T at one-loop from 
user-specified generic lagrangians (with approximate treatment of some of the loop correc- 
tions, see [9]) are currently available in micrOMEGAs. 

3 Indirect detection 

Indirect detection methods are many and varied. Here we focus on the following traditional 
methods: neutrinos from the Sun and the Earth, and gamma-rays, neutrinos, and charged 
cosmic rays (positrons, antiprotons and antideuterons) from annihilations in the galactic halo. 
There are also other indirect signals, like synchrotron emission, signals from cosmological 
halos (giving a diffuse flux), and indirect consequences of the presence of dark matter in 
stars, Chapter 29 [27], but we will not focus on them here. Most of the theory needed for 
this discussion is found in Chapters 24 [28] , 25 [29] and 26 [30] ; we use the notation in those 
chapters and elaborate on the formulae given there when needed. 

The main public tools available to calculate indirect rates are DarkSUSY [3j |4] [5] and 
micrOMEGAs [6] [Jj HJ [9] . In addition, there are also approximate simple formulae and param- 
eterizations available that can be used, but we will here focus on the numerical codes. 

3.1 Neutrinos from the Sun/Earth 

To calculate the neutrinos from the Sun/Earth, we need to calculate the capture rate of 
neutralinos in the Sun/Earth, we then need to solve the evolution equation for capture and 
annihilation in the Sun/Earth, let the neutralinos annihilate in the centre of the Sun/Earth 
to produce neutrinos and finally let the neutrinos propagate to the neutrino detector at Earth 
(taking interactions and oscillations into account). 

In Chapter 25 [29] . approximate formulae are given for the capture rate in the Sun. 
These formulae are good for quick calculations, but they include several approximations; 
with numerical codes, we can actually do better. DarkSUSY is currently the only public code 
that includes neutrino fluxes from the Sun/Earth and uses the full expressions in Ref. [51] . 
where the capture is integrated over the full Sun/Earth including capture on the 16 main 
elements for the Sun (and 11 for the Earth). In DarkSUSY, an arbitrary velocity distribution 
can be used if desired in place of the commonly assumed Maxwell-Boltzmann distribution. 
For example, the Earth does not capture WIMPs directly from the galactic halo, instead it 
captures from a distribution that has diffused into the solar system by gravitation interactions 
[32]. DarkSUSY uses a velocity distribution at the Earth based on numerical simulations that 
take this diffusion into account [33] . There are also indications from more recent numerical 
simulations of WIMP diffusion in the solar system [31] that heavier WIMPs will have a 
reduced capture rate in the Sun due to gravitational effects due to Jupiter. The DarkSUSY 
user can optionally include these effects. 

After capture, the evolution equation for the number density accumulated in the Sun/Earth 
is solved to give the annihilation rate today. Once the WIMPs have accumulated in the cen- 
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tre of the Sun/Earth, they annihilate and eventually produce neutrinos. In DarkSUSY, the 
annihilation and propagation of neutrinos is handled by a separate code, WimpSim [35j 136] . 
WimpSim takes care of annihilations to standard model particles in the central regions of the 
Sun/Earth with the help of Pythia [37]. Energy losses and stopping of particles in the dense 
environments at the center of Sun/Earth are also included. All flavours of neutrinos (and 
antineutrinos) are then propagated out of the Sun/Earth, taking oscillations and interactions 
(the latter only relevant for the Sun) into account. This is done in a full three-neutrino- flavour 
setup [35]. Once at the detector, the neutrinos are let to interact and produce charged lep- 
tons and hadronic showers. WimpSim has been run for a range of annihilation channels and 
masses from 10 GeV to 10 TeV, and the results have been summarized as yield tables that 
are read and interpolated by DarkSUSY. These results agree very well with a similar anal- 
ysis in Ref. [38], where parameterizations and downloadable data files are also given. For 
annihilation channels that are particle physics model dependent (like annihilation to Higgs 
bosons), the Higgs bosons are let to decay in flight in DarkSUSY and the resulting fluxes are 
calculated from their decay products (properly Lorentz boosted). 

The routines in DarkSUSY are also general enough to be easily adapted to other particle 
dark matter candidates, like Kaluza Klein dark matter. 

3.2 Charged cosmic rays 

The theory behind propagation of charged cosmic rays in the galaxy is presented in Chap- 
ter 26 [30]. We will use the notation in that chapter. In principle, what we have to do is to 
solve the master equation [30] Eq. (26.4)] with appropriate diffusion coefficients, energy loss 
terms, source terms and boundary conditions. Currently, micrOMEGAs includes the source 
spectra for arbitrary SUSY models, but does not include the spectra after propagation. This 
will be included in future versions though, using the results in Ref. [39 . In DarkSUSY, both 
the source spectra and the spectra after propagation in various propagation models are in- 
cluded. DarkSUSY implements axisymmetric propagation models and spherically symmetric 
(or at least axisymmetric) halo models. In DarkSUSY, diffusion is assumed to take place 
only in space (i.e. the term Kee in [201 Eq. (26.4)]) is assumed to be negligible). However, 
it also offers a full interface and integration with the leading cosmic ray propagation code 
GALPRDP [40] where more sophisticated propagation models can be used. For the halo den- 
sity, several preset profiles are available, the default being an NFW profile [41]. However, 
the user can supply her/his own halo profile; if it is given in the form of [30l Eq. (26.15)] it 
is particularly simple to do so. For solar modulation, DarkSUSY offers a standard spherical 
force-field approximation as explained in Chapter 26 [30]. 

For the source spectra (i.e. the spectra before propagation), DarkSUSY uses a similar 
setup as for the neutrino fluxes from the Sun/Earth given above. A large set of annihilation 
channels are simulated with Pythia [37 in vacuum, for a range of masses, and the yields 
of antiprotons, positrons, gamma-rays and neutrinos are stored as data tables. These tables 
are then read and interpolated by DarkSUSY at run-time. Higgs boson decays are included 
stepping down the decay chain. micrOMEGAs currently uses the same data files as DarkSUSY, 
but both codes are planning on using new updated simulations in future releases. 

For more details about the DarkSUSY implementation of antiproton fluxes, see Ref. [42]; for 
the positron fluxes, see Ref. [43]. The antideuteron fluxes are calculated from the antiproton 
fluxes with the method given in Ref. |44) . 

3.3 Gamma-rays and neutrinos 

Gamma-ray and neutrino spectra from annihilation in the galactic halo can be calculated 
with both DarkSUSY and micrOMEGAs. As mentioned above, the DarkSUSY spectra are 
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based upon Pythia simulations, which are then read in and interpolated by DarkSUSY. 
Currently, micrOMEGAs uses the same tables as DarkSUSY. DarkSUSY also includes internal 
bremsstrahlung photons [45] that can be very important in some parts of the parameter 
space (e.g. in the stau coannihilation region). DarkSUSY also includes the monochromatic 
gamma-ray lines from annihilation to 77 and Z7 [37] that occur at loop-level. Cur- 
rently, the DarkSUSY (and consequently then also the micrOMEGAs) neutrino spectra from 
halo annihilations do not include neutrino oscillations, but this will be addressed in future 
versions. 

Gamma-rays and neutrinos are not affected by propagation, and hence the flux can be 
written as 



dEdfl \ dE J VICT 29 crn^- 1 . 

xJ(ip) cm" 2 s" 1 sr" 1 , (20) 



where we have defined the dimensionless function 

'<*> = sk; - (osmsf)' im (21) 

with ip being the angle from the galactic centre direction to the direction of observation, 
and dN^/v/dE being the spectrum of gamma-rays or neutrinos. The line of sight inte- 
gral, Eq. (|2ip , can be calculated for any spherically symmetric profile in DarkSUSY (whereas 
micrOMEGAs currently implements an isothermal profile only). 



4 Exploring the parameter space 

One problem that arises when exploring a specific supersymmetric model setup (e.g. mSUGRA 
or a low-energy MSSM model) is how to scan the parameter space. The dimensions of this 
parameter space, i.e. the number of free independent parameters, can be large. For example, 
the general MSSM model has 124 parameters (MSSM- 124), 18 of which define the Stan- 
dard Model (SM). If one assumes CP-conservation, the number of parameters reduces to 63 
(MSSM-63), which is still a large number. Typically, one reduces the number of parameters 
still further with inspired theoretical insights (see Chapter 8 [48] ) . In Minimal Supergravity, 
for example, unification of coupling constants, of gaugino masses, and of scalar masses leaves 
only 23 parameters, i.e. SM+5, one of which is just a positive or negative sign. An interme- 
diate model often used in neutralino dark matter studies (MSSM-25) has 25 parameters, i.e. 
SM+7. 

One typically wants to find parameter values that are theoretically consistent, have a 
preferred relic density and are not already excluded by other searches (e.g. rare decays or other 
accelerator searches). The brute force method would be to scan over the parameter space 
with some kind of grid scan. One soon realizes that one can typically get a better efficiency 
in the scans (i.e. more points that pass the cuts, or a better sampling of different interesting 
regions in the parameter space) by scanning in the logarithm of the mass parameters instead 
of the mass parameters directly. For higher-dimensional parameter spaces, it is often also 
more advantageous to scan randomly instead of on a fixed grid. 

However, none of these methods are very effective in finding regions of parameter space 
that pass all the cuts. The relic density cut alone discards many models because of the high 
precision with which we know the relic density of dark matter today. Hence, more sophis- 
ticated methods have evolved that are efficient in generating points inside the interesting 
regions. Most of these use a Markov Chain Monte Carlo (MCMC) @9j W, EH EH E3J [54] to 
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Figure 2: A humorous scatter plot showing correlations that may arise with appropriate 
choices of priors while scanning parameter space: in this case, the priors reflect the anthropic 
principle. 



generate points according to a goal distribution specified by the researcher. The goal distri- 
bution can be a function without direct physical meaning (e.g. a Gaussian distribution that 
peaks in the desired region) or could have a statistical meaning (e.g. a likelihood function 
or a prior distribution in a Bayesian analysis). A recent public code to perform these tasks, 
and that is linked to DarkSUSY, micrOMEGAs and other codes, is SuperBayeS [53]. These 
advanced methods can be very effective in finding the interesting regions of the parameters 
space. However, when interpreting the distribution of points these methods produce, one 
should be very careful. The definition of "interesting" is different for different investigators, 
and the way points are generated always involves a prior in parameter space (even grid meth- 
ods can be said to have a prior, namely a series of Dirac delta functions at each grid point). 
One could go to the extreme of producing any kind of results by choosing appropriate priors. 
Fig. [21 for example, shows the correlation between the direct detection rate and the relic 
density jokingly obtained with priors that reflect the anthropic principle. When parameter 
space scans with priors are used to compute statistical inferences on data, e.g. likelihood con- 
tours [55] . one should keep in mind their rather severe dependence on the priors, especially 
the priors on very poorly known parameters for which there are little or no experimental data 
(supersymmetric masses, e.g.). This dependence arises from the use of Bayes theorem, which 
gives the probability of the model parameter given the data (the likelihood) in terms of the 
probability of the data given the model parameters (the assumed distribution of experimental 
and theoretical errors) as 

Prob(modelldata) = Prob(datalmodel) Prob ( model ) _ ( 2 2) 

Prob(data) 

While the probability of the data Prob(data) is just a normalization constant, the probability 
of the model parameters Prob (model) is the prior representing the degree of belief or the 
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Figure 3: Modules that are involved in the construction of a Monte-Carlo event generator 
for particle collider physics. 



relative preference the researcher has in specific values of the model parameters. When real 
experimental data on particle dark matter models are in, the dependence on the priors is 
expected to become less severe. 



5 Interface with collider and precision measurements 
codes 

Until a few years ago, one used constraints on the inferred amount of dark matter to delimit 
the parameter space of supersymmetry in order to narrow searches of supersymmetry at the 
colliders. Now with the improvement on the precision on the cosmological parameters one 
asks whether the LHC and LC can match the precision of the upcoming cosmology and as- 
trophysics experiments in, for example, reconstructing the relic density once supersymmetry 
is identified |56[ 154] . This may even bring a bonus in that one can test some cosmological 
and astrophysical assumptions, like indirectly "probing" the history of the early universe. 
For this programme to be feasible one needs to control the particle physics component with 
as much accuracy as possible. To be able to conduct a cohesive and self-consistent precision 
test of the origin of DM from the particle physics point of view one needs to calculate not 
only those dark matter cross sections but also the observables at the colliders that are pre- 
dicted for the same dark matter model. Ideally, therefore, one would like a common tool that 
performs this task. In many instances this also requires that one goes beyond calculations 
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Figure 4: A collection of codes for the calculation of observables for supersymmetric particles. 
The event generators and the tree-level matrix elements codes cover also the Standard Model 
and in some instances other models of New Physics. The arrows show how the output of one 
code is fed into another code, which can itself be a source for other codes. 



at tree-level. This is especially true in the case of supersymmetry where it is known that ra- 
diative corrections can be large. Some progress in this direction has also been made recently 
within the SloopS Collaboration [571 |5S1 159 ) 15U 1 [61 ] . 

The dark matter codes for supersymmetry such as DarkSUSY, IsaRED/RES and micrOMEGAs 
include some higher order effects. Moreover, because of the complexity of the MSSM which 
has a large number of parameters and a large array of predictions, these codes also rely on 
other more specific codes that predict various other observables in supersymmetry. This 
concerns for example codes for the calculation of the spectrum based on the renormalisation 
group equations (RGE) that predict the low energy physical masses from an input at the 
unification scale in some constrained model of supersymmetry breaking. Spectrum calcula- 
tors include Spheno [55], Softsusy [53], suspect [5J] and ISASUGRA (part of ISA JET) [TU] . 
These codes themselves may borrow from more specialised codes like those for the calcula- 
tion of the Higgs masses, such as Feynhiggs 55] . The codes for the mass spectra may also 
feed in stand-alone codes for the calculation of precision measurements, like the calculation of 
(g — 2)^ and B observables. Example of such codes or "flavour calculators" arc (SUSYbsg [66] . 
superiso [57]). 

To make contact with LHC and LC observables one of course needs matrix elements for 
cross sections and decays of the supersymmetric particles. Many multipurpose matrix element 
generators exist for supersymmetry. Among them, Amegic 68 , CalcHEP [16] . CompHEP [18] , 
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Grace-SUSY [69], Omega [70], Madgraph [71]. Multipurpose matrix element generators can 
return results for any cross section or decay of a supersymmetric particle at the tree level. 
More dedicated and specialised codes in this category (cross sections and decays) usually 
improve by going beyond tree-level, among them PR0SPIN0 [72] for the production of su- 
perparticles at a hadronic collider, HDECAY [73] for the decay of the Higgs and SDECAY [73] 
for the decay of other sparticles. Automatic codes for generic one-loop cross sections and 
decays with supersymmetry have also been completed recently as concerns the electroweak 
corrections and some classes of QCD corrections: Grace-susy-lloop 75] and SloopS [59]l60]. 

For simulations at the colliders one still needs to incorporate the matrix elements for the 
production (the hard process) and the decays (of the unstable superparticle resonances) into 
fully fledged Monte-Carlo generators. The latter include (i) parton shower (radiation), (ii) 
multiple interaction and beam remnant in hadronic machines and (iii) hadronisation. The 
main Monte-Carlo event generators are currently Herwig++ [76], ISAJET [10], Pythia [37] 
and SHERPA [77]. Fig. [3j shows the ingredients that go into a Monte-Carlo event generator. 
The mass spectrum module is what defines the model here, its content is therefore also 
encoded in the dark matter codes. One can use these generators to simulate the signatures 
of a particular model at the colliders and combine these findings with the manifestation of 
the same model in dark matter searches (direct and/or indirect), the prediction of the relic 
density in a particular cosmological model. One can constrain or reconstruct the model even 
further by taking into account observables from indirect precision measurements encoded in 
the flavour calculators. Codes ("Fitters") that perform these fits or constraints have been 
written specifically with supersymmetry in mind; we can mention Fittino [78], SFitter [79] 
and SuperBayeS [53] . 

As we have seen, there is a very large variety of codes that cover different aspects of the 
phenomenology of supersymmetry. A recent compendium of these codes can be found in [80] . 
Because of the large number of parameters in a general supersymmetric model and because 
many modules are fed into other modules, it is best to avoid errors as much as possible when 
passing parameters from one code to another. Some of these errors can be as trivial as a 
problem of sign convention. The SUSY Les Houches Accord (SLHA) jT3] allows an easy 
parsing for the MSSM parameters. This accord has been extended [14] to deal with more 
general supersymmetric models, like the NMSSM for which a version of micrOMEGAs exists 
that uses or can be used with the NMHDECAY code [81] . or with the inclusion of CP violation 
via CPSuperH [82] . 

Fig. 0] shows how the different codes for the calculation of supersymmetric observables 
both at the colliders and in dark matter searches or the evaluation of the relic density are 
interrelated. The calculation of the matrix elements needed for these codes requires first of 
all reading the Feynman rules. This in itself is a titanic endeavour because of the complexity 
of the MSSM and its extensions. For the dark matter codes, where a very large number of 
processes are involved, especially in the calculation of the relic density, practically the whole 
set of rules is called for. This is even more so for one-loop calculations. Special tools now 
exist to achieve this task whereby the model file (containing the Feynman rules) is generated 
automatically from just coding the Lagrangian in a manner as close to the calculation by 
hand. This was first done more than 10 years ago by LanHEP [17] based on C for easy 
interface with CompHEP. The implementation here is very similar to the canonical coordinate 
representation. Use of multiplets and the superpotential is built-in to minimize human error. 
The ghost Lagrangian is derived directly from the BRST transformations. Very recently 
FeynRules |83] based on Mathematica enables to perform the same task and can output to 
different matrix element generators. MicrOMEGAs rapid development was in large part made 
possible because of the extensive automation based on LanHEP and CalcHEP/ CompHEP. LanHEP 
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has also been greatly extended to automatically implement a model file for calculations at one- 
loop. It is now possible to shift fields and parameters and thus generate counterterms most 
efficiently. LanHEP has been successfully interfaced with the highly efficient and automatised 
one-loop packages based on Feynarts/FeynCalc/LoopTools 84, 85, 86] [87] . The SloopS 
package is the combination of LanHEP and FeynArts, after a fully consistent and complete 
renormalisation of the MSSM has been completed [57 , 58, 59, 60 . SloopS has been developed 
from the outset such that it is applicable to both high energy collider observables and for 
processes occurring at very low velocity such as is the case with dark matter particles. This 
will bring the cross breeding between dark matter and the collider predictions to a new level 
of precision, at least as far as the particle physics component is concerned. 

Combining codes to be able to conduct global analyses for dark matter searches and 
the determination of its microscopic properties are witnessing an intense activity. If future 
colliders discover SUSY particles and probe their properties, one could predict the dark 
matter density and would constrain cosmology with the help of precision data provided by 
WMAP and PLANCK. It would be highly exciting if the precision reconstruction of the 
relic density from observables at the colliders does not match PLANCK 's determination. 
This would mean that the post-inflation era is most probably not radiation dominated (see 
Chapter 7 [11] for a discussion of alternative cosmologies before Big Bang nucleosynthesis and 
their effect on particle dark matter). The same collider data on the microscopic properties 
of DM, when put against a combination of data from direct /indirect detection, can also 
give strong constraints on the astrophysical properties of DM such as its distribution and 
clustering. These properties reveal much about galaxy formation 56, 54j. For this program 
to be carried through successfully, tools developed in the cosmology /astrophysics community 
and tools developed within the particle physics collider community (and their interfaces) are 
essential. 
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